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Abstract 

We present a method to compute, by numerical simulations of lattice QCD, the 
inclusive semileptonic differential decay rates of heavy hadrons and the structure func- 
tions which occur in deep inelastic scattering. The method is based on first prin- 
ciples and does not require any model assumption. It allows the prediction of the 
differential rate in B semileptonic decays for values of the recoiling hadronic mass 
W ~ y/MghQCDi which is in the relevant region to extract \V u b\ from the end-point 
of the lepton spectrum in inclusive decays. 



1 Introduction 



In this paper we propose a method to compute, on the lattice, the shape function 
f(k+) which enters the calculation of the inclusive differential semileptonic decay 
rate of heavy hadrons. The knowledge of f(k+) is a fundamental ingredient for the 
extraction of \ V u b\ from the end-point of the lepton spectrum. The same function also 
enters the calculation of the photon spectrum in radiative B decays. Furthermore, 
the same method can be applied to the calculation of the structure functions of Deep 
Inelastic Scattering (DIS). Our approach does not require any model assumption and is 
based on standard techniques which are used to compute the hadronic matrix elements 
in lattice QCD. In particular, we show that the use of the Euclidean space-time, which 
is unavoidable in lattice calculations, is not an obstacle to the calculation of f(k+) in 
the deep inelastic region. 

Our main result is that, from the study of suitable combinations of lattice Green 
functions, we can obtain the quantity 



e -VQ 2 t [M B _ 

G{t,Q) = —= r f(k + )e~ k +\ (1) 

2VQ 2 Jo 

where Mb is the mass of the decaying hadron (the B meson in the example considered 
in this paper), t is a time distance and Q a spatial momentum that we can inject in 
the Green functions computed in the numerical simulations. By varying t and Q, we 
can unfold the integral and extract f(k + ). 

The function f(k+) in eq. (|l|) differs from the usual shape function f(k+) introduced 
in refs. §], as will be explained in the following. In terms of f{k + ), the differential 
semileptonic decay rate can be written as 

dT - [ MB ^ f/u ^ dTpM < * fiM 

m = I HeT b 

q2 i-Mb r 

= I^I'y^s E t J dk + tt k +) @ « ~ 2E ^ [ 3m &* 2 " 4 < E 4 ' ( 2 ) 

where m£ = Mb — k+. The advantage of using f(k + ) is that any reference to unphysical 
quantities, such as A or the quark mass mi,, disappears in the expression above and 
the differential rate is written in terms of hadronic quantities only. 

The paper is organised as follows. In section Q we review the main formulae for 
the total and differential semileptonic decay rates as they can be derived using the 
Operator Product Expansion (OPE); f(k + ) is introduced in this section. In sec. || we 
present the basic formalism needed to explain our idea. We discuss, as a prototype, 
the case of the shape function which enters the differential decay rate, and then extend 
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the method to the computation of the structure functions. We express these quantities 
in terms of suitable T-products of local operators and discuss their continuation to 
the Euclidean space-time. In sec. we describe the implementation of the method 
in lattice calculations and discuss the systematic errors which may arise in actual 
numerical simulations, due to the present limitations in computer resources. 



2 Inclusive decay rates 



This section contains a summary of the main results obtained in the literature for 
semileptonic and radiative inclusive decays. We present below the main formulae for 
the total and differential rates which can be derived using the OPE This will 

allow us to introduce the definition of the shape function and of its moments. 

The idea that inclusive decay rates of hadrons containing heavy quarks can be com- 
puted in the parton model is quite old and was used, for example, in ref. ||] to predict 
the charmed-hadron lifetimes. To account for bound state effects, the partonic calcu- 
lation was subsequently improved by the introduction of a phenomenological model, 



called the "spectator model" [|H], 11], which has been and continues to be extensively 
used to extract \V c b\ and \V u b\ from inclusive semileptonic decays. A noble theoretical 
framework for the spectator model was then provided by the the Wilson OPE applied 
to the inclusive decays of heavy hadrons 1|. The expansion parameter is not neces- 
sarily the inverse heavy-quark mass I/uiq, rather the inverse of the energy release 1/W 
of the process at hand. In a large region of the available phase space W is of the order of 
niQ. In this case, under the hypothesis of quark-hadron duality, the operator-product 
expansion is expected to give accurate predictions for the decay widths, expressed in 
terms of few non-perturbative parameters. In particular, at lowest order in l/m,Q, the 
expression of the decay widths derived with the OPE coincides with that obtained in 
the parton model calculation. 

Let us recall how this works for the inclusive semileptonic processes B — > X CtU £uji. 
Using the OPE, the total semileptonic decay rate is given by Q, |5| 



T(B — > X CtU £i>() = |V( CjU )b 



1927T 3 



(3) 



where we have neglected QCD radiative corrections and phase space effects due to the 
mass of the final charged lepton. The dots stand for higher-order perturbative and/or 
power corrections. We have omitted to write explicitly the terms of 0(l/m|) because 
they are expected to give negligible contributions to r(.B — > X cu £vi^. Co ) i(^c,w) are 
known phase-space factors which depend on the ratios of the final-quark masses to nib, 
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x c ,u = ^cu/ 17 ^ @> |I- I n e 1- ©) the term which coincides with the parton model 
result 

T PM (B -> X C)U ^) = |y (CjU)fe | 2 C (x c , M ) (4) 

corresponds to the insertion of the leading, dimension-three operator bb appearing in 
the OPE. In the heavy quark effective theory (HQET) at fixed velocity v, this is a 
conserved operator (66 — > h v h v + . . .) satisfying the following normalisation condition 

(B(v)\h v h v \B(v)) = 1. (5) 

The hadron state \B(v)} is normalised to v° instead of the usual relativistic normal- 
isation of 2E because this is more convenient for the heavy quark expansion. The 
corrections of 0(1/ mi) vanish since the only possible dimension-four operator, bilpb, 
can be reduced to 66 by using the equations of motion. The kinetic energy of the 
heavy quark in the B meson, Ai, and the chromomagnetic moment of the heavy quark, 
A2, correspond to matrix elements of local, dimension- five operators appearing in the 
OPE pj. Written in terms of the fields of the HQET, these operators are given by 

Ai = (B(v)\h v (iDfh v \B(v)) 

A 2 = (B(v)\h v a^iDHD v h v \B(v)) . (6) 

Thus, up to 0(l/mf) the problem is reduced to the calculation of the matrix elements 
of the two local operators appearing in eq. (0). For these matrix elements several 
theoretical estimates exist 0. A recent review of these estimates can be found, for 
example, in ref. [p^ j. 

Unfortunately, in order to determine \V u b\ from the lepton spectrum, one has 
to study the differential distribution and perform cuts to suppress decays involving 
charmed particles in the final state. For the differential distribution, the relevant scale 
of the OPE is the squared momentum of the hadronic system recoiling against the lep- 
ton pair, Qt, = rribV—q, where q is the momentum of the lepton system. By approaching 
the end-point of the electron spectrum W = J Q\ becomes small. 

In the region where v ■ Qb is of 0(Mb) and W ~ \/M^Kqcd we can still use the 
parton model by introducing a "shape function" , analogous to the distribution function 
of the spectator model |Ti"|| , corresponding to a modified OPE. This is analogous to deep 
inelastic scattering, where one can introduce non-perturbative distribution functions 

1 We are not concerned here with the precise definition of Ai ; 2 when QCD corrections are taken into 
account. 
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for the light-cone component of the quark momenta. In terms of the shape function 
f(k+) the differential distribution is given by 

dT f A „ .„ s dT 



dk + f(k + ) ^L( m *,E e ) 



dEi J-m b ' dE, 
G 2 - ^ 



= \ V ub\ 2 j^ El j ^ dk + f(k+) Q(m* b - 2E e ) [3m* b 2 - 4m* b Et] , (7) 

where m b = m b + k + and 

f(k+) = (B(v)\ h v S(k+ - iD+)h v \B(v)) . (8) 

D + = n ■ D, with = Q b /(v ■ Qb)- Note that in the region where W ~ ^JMb^-qcd^ 
n 2 ~ Aqcd/Mb ~ 0. 

Finally, if we consider the extreme region where W ~ Aqcd, the partonic picture 
breaks down. In this region the rate is dominated by few single states or resonances 
and the appropriate description can only be obtained by predicting the rate as a sum 
over the lowest exclusive channels. 

At the lowest order in l/m b , the integration of dT/dEi over Eg in eq. (0) gives the 
total rate. It also gives some of the 0(l/m 2 ) corrections appearing in eq. @. For the 
following discussion, it is convenient to introduce the moments of f(k + ), defined as 

M n = [ A dk + k n + f{k + ) . (9) 
J-m b 

The first few moments are known: the leading contribution, given by Mq = 1, is the 
normalisation of /(&+) and is fixed by the normalisation of the conserved scalar den- 
sity operator of eq. (||); the second moment M.\, related to the correction of 0(l/rrih), 
vanishes because there are no operators contributing at this order; the correction of 
0(l/m 2 ) is given by M.2 = — Ai/3. The contribution of the heavy-quark chromomag- 
netic moment is not included in the shape function. 

In the case of the integrated rate, the contributions coming from higher moments 
of the shape functions are related to operators of dimension larger than five. For this 
reason, they are suppressed by higher powers in l/mb an d can be safely ignored. For 
the differential distribution of eq. (0) near the end-point, instead, all moments 
become of the same order in 1 /m;, and cannot be neglected 0. The knowledge of the 
full shape function /(/c+), and not only of the first few moments, is then needed in this 
region. This is analogous to what happens for the structure functions in DIS processes 
in the region where the Bjorken variable x — ► 1. 



A more extensive discussion of this point will be presented in section ^. 
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So far, several models for the shape function have been proposed 0, 11, 13]. In 
particular, one can show that the old ACCMM model is equivalent, up to a process- 
dependent redefinition of the heavy quark mass Q, to a specific choice of the shape 
function. It remains true, however, that the extraction of \V u f,\ from the experimental 
measurements of the end-point of the lepton spectrum using these models is "model 
dependent" and consequently affected by a systematic error which is difficult to esti- 
mate. 



3 The shape function f(k+) 



In this section, we show that it is possible to determine the full shape function, and 
not only a few moments of it (as for example Ai |l4j|), by numerical simulations on the 
lattice. The shape function can be obtained from the time and momentum dependence 
of suitable Euclidean Green functions of the same type as those which are currently used 
on the lattice to compute hadronic matrix elements of local operators. We show that, 
in the inelastic region where the knowledge of the full shape function is necessary, this 
can be done in spite of the fact that lattice simulations are performed in the Euclidean 
space-time. This is new since, as discussed below and in ref. |15|], the use of the 
Euclidean space-time, instead of the Minkowsky one, seems to prevent the possibility 
of computing inclusive quantities on the lattice 0. We also show that, by using the 
same approach, it is possible to compute the structure functions of DIS, and not only 



few moments of it, as done so far [17]— yi 

A further advantage of our proposal is the following. Higher moments of the shape 
function (and those of the structure functions) are related to matrix elements of local 
operators of higher dimensions which are, in general, afflicted by power divergences 
when a hard cutoff is used (or renormalon ambiguities in dimensional regularisations). 
This makes the definition of the renormalised operators, which have finite matrix ele- 
ments when the cutoff is removed, very problematic [jnj. With our method, instead, 
no renormalisation is needed as we get directly the shape function from the lattice 
correlation functions, as explained in sec. ^. 

We now give all the details of the derivation of our method. 



3 For exclusive processes involving more than one particle in the final state, the problem is usually referred 



to as the Maiani- Testa no-go theorem 16 . 
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The spin- averaged, differential semileptonic decay rate is given by 



dq 2 dE f dE Ue 1 (c ' M)b| 2vr3 



IT, + W 2 ( 2£^ - X -q 2 ) + IT <f it) - ) 



(10) 

with the Lorentz-invariant functions Wj defined by the expansion of the hadronic tensor 
= (27rfJ2^(PB-q-Px)(B(v)\J^\X)(X\r\B(v)), (11) 

X 

which can be written as 

W p.v = -gi^Wi + v^v v W2-id wc ^q a VfiW3 

+ gY^ + fgV + gV)^. (12) 

In eq. (|Tl|) , the weak current is defined as J M = q(x)^b(x) , where q{x) represents the 
field of the light final quark and 7^ = 7^(1 — 7s)/2. Using the optical theorem, the 
form factors W{ are given by 

Wi = --hnTi , (13) 

7T 

where the T are defined by the forward matrix element of the T-product of the two 
weak currents 

= -i J d 4 x e~ iq - x (B(v)\T (> \x)J v {Q)) \B(v)} 
= -gTTx + v»v»T 2 - ie^ af3 q a vpT 3 

+ q»(fTi + {<fv v +q"v»)T B . (14) 

Thus the problem is reduced to the calculation of the forward matrix element of the 
T-product defined above. 

One could think that, with sufficient computer resources to work with very small 
values lattice spacing a such that 

Aqcd <C m b < - , (15) 

it would be possible to compute directly W^ v on the lattice using the T-product of 
eq. (|j~4|). This is in general illusory. The problem arises from the fact that numerical 
simulations are necessarily performed in the Euclidean space-time Let us illustrate 
where the problem comes from. 

In the Minkowsky case, for t > 0, one has 



W^ u = -Im 

7T 



i / d 4 x e~ iq - x (B(v)\J^ \x) J u (0)\B(v)) 
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— Im 

7T 



W dt e «&-4°-**)* (27r) 3 5 3 fe _ f_ px) x 

(B(v)\J^(0)\X)(X\r(0)\B(v)) 

= Y,( 2 ^f 6 HPB-q-px)(m\J tl H0)\X)(X\r(0)\B(v)) . (16) 
x 

For t < 0, the cut corresponds to an intermediate state with two 6-quarks and a c, which 
is not the process we are interested in. The 5 function related to energy conservation 
is given by the integral over the time of an oscillating exponential. In the Euclidean 
case, instead, eq. ( jl6| ) becomes 

W>* W ~J2[ dt e^B-tP-Ex* (2n) 3 S 3 (p B -q-px)(B{v)\J^^{0)\X)(X\r(0)\B{v)) . 
x J 

(17) 

Thus, at large time distances, the integral over t is dominated by the states with the 
smallest energy Ex, and not by those which satisfy the energy conservation condition 
Ex = P% ~ 1 ■ This is similar to what happens in the case of exclusive decays which 
was discussed in ref. (l(|. The same problem is also encountered in other analytic 
approaches which study inclusive cross-sections or decay rates in the Euclidean [20|. 

We now show that, in the deep inelastic limit, the analytic continuation to the 
Euclidean space-time of the T-product of the currents in eq. fll4| ) becomes very simple 
and that it is possible to extract the wave function by studying suitable Euclidean 
Green functions. 

T^ u can be written as 

T*" = -ijd*x e~ iq - x {B{v)\b{x)^ L S{x - 0h£b(0)|5(u)> 

= -i J d*x j(PB-i)-* {B{v)\b v {x)^ L S{x - 0)-/Zb v (0)\B(v)) 

-i I d 4 x {B(v)\b v (x) 7 £S Q (x-0hZb v (0)\B(v)) , (18) 



where S(x — 0) is the final light-quark propagator; b(x) = e~ ipB ' x b v {x) and Sq(x — Q) = 
e l Q x S[x — 0), with Q = pb — Q = Mbv — q. Our definition of the heavy quark field 
b v (x), written in terms of the physical momentum of the meson ps instead of the quark 
momentum p&, is legitimate since the shape function is only defined at the lowest order 
of the heavy quark expansion. The reason for this choice will be explained below. It 
is possible to expand Sq(x — 0) as follows (for clarity we neglect the mass of the light 
quark) 

Sq(x-O) = 



+ 4p + ieJ x _ \Q 2 + 2iQ ■ D - D 2 - l/2a^G^ + ie ) x _ Q 
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Q 2 + 2iQ ■ D + ie 



(19) 

x-0 

In eq. ( |l9| ) we have kept only the leading terms of the expansion in powers of 1/ \J Q 2 = 
1/W ~ 1/Mb, and those which become leading near the end-point of the lepton 
spectrum, i.e. where W ~ MbAqqd - f|- Using eq. (19) one finds 

1 / \ r M B f(b , ) 



1 / \ r M B 

= - (QV + Q v v» - gTQ ■ v + ie^Q a v p ) J dk A 



Q 2 -2v Qk + + ie ' 
(20) 



where the shape function /(&+) is defined through the relations 
(-l) n (B(v)\b v j u (W^) . . . (iD^)b v \B(v)} = M n v u v^ . . . v^+BnS^v^ 2 . . . v" n +. 

The moments Ai n in the equation above are given by 



(21) 



rM B 

M n = / dk + k n + f(k + ) (22) 
Jo 

With our choice of Q, written in terms of the 5-meson momentum, and up to a trivial 
change in sign, the k + of eq. (0) is changed into — k + + A, leading to eq. (§). Note that 
with our definition of k+, ml = Mb — k + , which is in our opinion more physical: the 
differential distribution is now expressed in terms of hadronic quantities only, without 
any reference to unphysical quantities such as the quark mass or A. 

The moments A4 n ~ ^qcd have been defined in eq. (^) and give rise to terms 
proportional to M n (v ■ Q/Q 2 ) n ~ (AqcdMb/W 2 ) 11 . Thus, their contribution to the 
rate is suppressed as (Aqcd/Mb) 11 when W ~ Mb, whereas it becomes of 0(1) in the 



region where W ~ y/ ' MbAqcd- The contributions proportional to B n are subleading 
in 1/Mb with respect to the corresponding A4 n in all the physical region of interest, 



i.e. for all values of W, including W ~ \/ MbAqcd (but not in the elastic region where 
W ~ JXqcd)- Note that in order to derive eq. (f20|), we have never used the HQET, but 
only the OPE, by separating the large frequency modes (~ W) from the low energy 
ones (~ Aqcd) and expanding in powers of Aqcd/W. For simplicity, without loss of 
generality, we consider in the following the case of a B meson at rest, namely v = (1,0). 
We now consider, for t > 0, the Fourier transform of T^ v defined as 

T^{t,Q) = J f^£ e -*3ot r /tt' _ (23) 
In terms of the T-product of the currents, we have 

T^(t,Q) = -if d 3 x (B(p B = 0)\J^(x,t)J»(0)\B(p B = 0)) 

= - ie -iM B t J d 3 x e -iQ-x (p($ B = 0)\jrt(x,t)J u (0)\B(p B = 0)),(24) 
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where Jtf(x) = q(x)^b v (x). The factor e ~* A/ s* appearing in eq. ( |24| ) cancels the 
corresponding term in the three-dimensional Fourier transform of the correlator, so 
that the r.h.s. of this equation goes as e - *^ *, which is the expected behaviour for a 
hadronic system with energy Qq. 

For t > 0, by closing the contour of the integration over Qq below the real axis, it 
is straightforward to find 

j rMg . _ _ . 

T^(t, Q) = -- J dk+ f(k+) (Q»6 u0 + Q u 5^° - g^Q° - ie°^ a Q a ) x 



2\IQ 2 

where Q = (Qq , Q), with 



(25) 



Q+ = k + + ^Jk\ + Q 2 ~k + + VQ 2 • (26) 

Using eq. (^3|) it is very easy to make the analytic continuation to the Euclidean 
space-time 

W^(t,Q) = --lmT^(t,Q) (27) 

7T 



= — / dk + f(k+)- = x (Q»5 u0 + Q"5^° - g^Q° - ie°^ a Q c 

By a suitable choice of the Lorentz components fi and v of the currents and of the 
spatial momentum Q, one can isolate either an integral of the form given in eq. (jl]) or 
the following one 



G(t, Q) = / dk + f(k + ) e- fc + 4 Q+ . (28) 

"2 Jo 



2\jQ 

By studying the time dependence of G(t, Q) at several values of Q, we can unfold both 
the integral above and the one in eq. (|]) and extract the shape function. 

One may be surprised that the analytic continuation to the Euclidean space-time 
is so simple, in spite of the argument made at the beginning of the section. Indeed, 
also in the Minkowsky case, in order to apply the OPE, one has first to expand the 
T-product of the currents for large space-like momenta, which is equivalent to work in 
the Euclidean, and then continue the resulting expression to the kinematic region of 
interest. If this is possible, which implicitly corresponds to the assumption of quark- 
hadron duality, then the expansion has the very simple form of eqs. (BO) and (E~~ 



9 



In these equations, the singularities are those of a free particle propagator, the Wick 
rotation of which is straightforward and leads to the result given in eq. (|27|) . The con- 
tamination of intermediate states with small invariant masses, which would dominate 
at large time distances, is suppressed by injecting in the correlation functions a large 
spatial momentum Q, with v/Q^ 3> Aqcd- 

We now show that the same method can be used to compute the structure functions 
of deep inelastic scattering. The starting formula is, as before, the T-product of two 
currents. We have 

T» v = -% J (fx e- iq ' x {M\T (> t(ar) j"(o)) \M) , (29) 

where \N) represents a generic hadronic state and q is the momentum of the external 
vector boson with space-like momentum q 2 < 0. In this case, we have only q as large 
momentum in the game. Thus we may write 

= -if (fx e~ iq - x {M\^{x)T^S{x - 0)PV(0)|A/) 

= -i J (fx (N\4){x)T^S q {x-V))T u i>{0)\N) , (30) 

where = 7^ or = 7^ for electromagnetic or weak charged currents respectively 
and S q (x — 0) = e~ iq ' x S(x — 0). By expanding S g (x — 0) in powers of — q 2 , as done before 
for the inclusive decay rate, we may define a distribution function also in this case. 
Thus the same formulae apply to the case in which the large momentum is space-like, 
as in DIS, or time-like, as in inclusive B decays. 



4 Implementation of the method in numerical 
simulations 

In this section, we briefly explain the method to extract f{k + ) from the Euclidean 
lattice correlation functions and discuss the feasibility of our approach. Without loss 
of generality, we work with a B meson at rest. The extension of our method to a B 
meson with an arbitrary velocity v is straightforward. 
W^ u (t, Q) can be readily obtained from the ratio 

W^ v (t Q) 

W^{t,Q)= lim W ' g~Mgt ^ (31) 

tf,ti->00 Otf,ti 

where 

W t %(t,Q) = lfd 3 x e- i ^(0|4 B=0 (t / )jt(^ t )j iy ( )$ p -. B=0 (-t i )|0), (32) 
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and 

S tf,U = <0|*L=o(*/)*?B=o(-*<)|0) • (33) 
$p- B (t) is the B interpolating field with definite spatial momentum p B 

$p B (t) = J d 3 xe- ipB - s ^ B (x,t) . (34) 

It is very easy to demonstrate eq. ([H]). In the Euclidean, using the transfer matrix 
formalism, we have 

$ B (x,t) = e At <f> B (x)e- At , (35) 

so that the correlation functions have an exponential dependence on the energy of the 
external states. This implies that, in the limit U,tf — > do, the lightest boson state, 
corresponding to a B meson, dominates the correlation functions (^) and (|33|), since 
all higher-energy states are exponentially suppressed. In this limit, with t > 0, we then 
obtain 

W% ti (t,Q) - 1(0|^(0)|B(pb = 0)) x [/ d 3 x e'^(B\4(x,t)J u (0)\B) x 

(B(p B = 0)|*fl(0)|0)e-^(** + ^ , (36) 

and 

St f , U - <0|$i(0)|5(p B = 0))(B(p B = 0)|$ B (0)|0)e-^^+*/) . (37) 

Thus we have shown that for U,tf — * oo the ratio in eq. ( |3l| ) directly gives the required 
quantity and that the knowledge of the coupling of the interpolating field $b to the 
physical meson state is not necessary. 

We now discuss the feasibility of the method in actual numerical simulations. If 
one uses the complete light quark propagator, a systematic effect in the extraction of 
f(k+) in the relevant region, k + ~ A, is induced by contributions of states with a small 
invariant mass, which dominate the Euclidean correlation functions when t — > oo, see 
eq. (|17|). By injecting a large momentum Q in the correlation functions, the difference 
between the energy Ex of these states and the energy Q$ of the states we are interested 
in becomes Qq — Ex ~ Qq — \J Q 2 ~ A. Moreover the contribution of these states is 
suppressed at least as A 2 /M B . In order to compute f(k + ), the condition 

T = ^e Kt « 1 (38) 
Ml 

must then be satisfied. Thus eq. (|3S|) gives a condition on the maximum time distance 
which can be used to extract f{k + ). We can formulate the condition on the maximum 
allowed time, in units of the lattice spacing, as follows 

2V t = -<^M#)- (39) 
a Aa A. A 
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By working with an arbitrarily small lattice spacing a, the condition (^) can always 
be satisfied for large values of JV* and f(k + ) can then be computed with negligible 
uncertainty. In practice, the value of a is limited by computer resources and we have 
to worry whether we have a sufficient number of points in time to unfold the shape 
function. To give an example, with a lattice spacing a -1 = 6 GeV, a value which will 
probably be reached with the 100-gigaflops/teraflops machines already available or in 
construction, a heavy meson mass M ~ 2 GeV and A = 0.3 GeV, one finds Nt ~ 75, 
which is certainly large enough. Note that f(k + ) cannot be studied for k + S> A because, 
in practice, the correlation function is dominated by subleading terms in 1/Mb at all 
accessible time distances. This region, however, can be studied perturbatively. 

In order to avoid the limitations on Nt, a possibility is, of course, to perform the 
numerical simulations using the approximate propagator of eq. (jl^) (we omit the $ in 
the numerator because this gives rise to trivial kinematic factors) 

^ Q = Q 2 + 2iQ ■ D + ie = (tTq) Q 2 /v ■ Q + 2iD + + ie ' ^ 
It is straightforward to show that Sq(x) can be written as 

e i(vQ)x+/2 

Sq(x) = — — q — S LEET (x) , (41) 

where x + = n ■ x and S LEET {x) is the light-cone propagator of the Large Energy 
Effective Theory (LEET) Q , which satisfies the equation 

2iD + S LEET {x) = 5\x) . (42) 

Thus, the extraction of f(k+) from Sq(x) is equivalent to the use of the LEET. In 
our case, we are allowed to use this approximation since it has recently been shown 



that the LEET is applicable to inclusive processes [ 24 ] , in spite of the difficulties that 



it may have for exclusive decays [25]. Note that the calculation of the physical shape 
function using the LEET propagator, which is more singular than the propagator of the 
full theory, requires a further logarithmic renormalisation of W^it, Q) |^], which can 
be computed in lattice perturbation theory. The ultraviolet divergences of the LEET 
correspond perturbatively to infrared divergences in the full theory |2^]. In the latter 
case the infrared divergences are automatically regularised by the non-perturbative 
contributions in the physical matrix elements of the T-product of the two currents and 
no renormalisation is required. 

It is not clear to us whether the use of S LEET (x) will be convenient in practice, since 
this propagator is much more singular than the full one. For this reason we expect 
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that the correlation functions computed in numerical simulations using S LEET (x) will 
be affected by larger statistical fluctuations. 

We stress that in all the formulae derived in this paper, we never used the fields 
of the HQET and, indeed, the same formalism also applies to the calculation of the 
structure functions in DIS. The reason is that the shape function is defined from the 
OPE in powers of 1/W (at lowest order in 1/toj,) and one has only to worry about those 
terms which become of 0(1) at the end-point of the lepton spectrum. Thus, the use of 
the fields of the HQET would only affect terms of higher order in 1/ mj on which we do 
not have control anyway. In the calculation of the relevant correlation functions, one 
may use the heavy-quark propagator of the lattice HQET for the b quark. We recall, 
however, that calculations done using the lattice HQET are afflicted by considerable 



difficulties, which are absent in the full theory [26, 27]. Thus, we believe that the best 
strategy is to compute the shape function in the full theory, both for the light and the 
heavy quarks, at several values of the heavy quark mass, and extrapolate the results to 
the B case. This strategy has been already successful in the calculation of the heavy 
meson decay constants, of the form factors for exclusive semileptonic decays, and of 
the B°-B° mixing ^-parameters pq, 27 1. 



5 Conclusions 

In this paper, we have shown that it is possible to compute, by numerical simulations 
of lattice QCD, the shape function which describes the lepton spectrum in semileptonic 
B decays and the photon spectrum in radiative B decays. The same approach can be 
used for the calculation of the structure functions in DIS. The method is based on 
first principles and does not require any model assumption. Moreover it avoids the 
calculation of the matrix elements of higher dimensional operators, which are plagued 



by power divergences and the renormalisation of which is very difficult to achieve [17]. 
Indeed, no renormalisation is needed and the shape function can be extracted directly 
from suitable ratios of lattice correlation functions. We have also proposed a redefini- 
tion of the shape function which avoids any reference to A or to the b-quark mass, and 
allows us to write the differential rate in terms of hadronic quantities only. 
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